Self-Replicating Three- Dimensional Vortices in Neutrally-Stable Stratified Rotating 

Shear Flows 
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A previously unknown instability creates space-filling lattices of 3D vortices in linearly-stable, 
rotating, stratified shear flows. The instability starts from an easily-excited critical layer. The 
layer intensifies by drawing energy from the background shear and rolls-up into vortices that excite 
new critical layers and vortices. The vortices self-similar ly replicate to create lattices of turbulent 
vortices. The vortices persist for all time. This self-replication occurs in stratified Couette flows 
and in the dead zones of protoplanetary disks where it can de-stabilize Keplerian flows. 
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Introduction. For a protostar to accrete gas from its 
surrounding protoplanetary disk (PPD) and grow into a 
star, the PPD must have an instability capable of trans- 
porting angular momentum radially outward [1]. This 
requirement has led to an intense effort to find insta- 
bilities in PPDs and other rotating flows in which the 
angular momentum satisfies Rayleigh's criterion for cen- 
trifugal stability, i.e., the absolute value of the angular 
momentum increases with increasing radius [H. Numer- 
ical studies Jsl, Q of ideal gases in PPDs and experimen- 
tal studies [5| of constant density fluids in rotating flows 
where the velocity obeys Rayleigh's criterion show that 
the flows are both linearly and finite- amplitude stable 
to hydrodynamic perturbations. In regions of a PPD 
where the gas is sufficiently ionized to couple to mag- 
netic fields, the magneto-rotational instability (MRI) [l| 
can operate. However, large regions of PPDs, known as 
dead zones^ are too cool to ionize sufficiently and have 
MRI. Strato-rotational instabilities [6] and Rossby wave 
instabilities could de-stabilize a PPD, but they re- 
quire unrealistic impermeable radial boundaries on the 
PPD or a large continually-forced density perturbation, 
respectively. Thus, star formation remains problematic. 

Here we report a new type of finite-amplitude insta- 
bility that occurs in neutrally-stable, shearing flows with 
stably-stratified densities (as in a PPD) and with veloc- 
ities that would satisfy Rayleigh's stability criterion if 
the densities were constant - as assumed in Rayleigh's 
analysis. These flows include plane and circular Cou- 
ette flows, which have been used to model PPDs. Our 
new instability allows a small-size perturbation to form 
large vortices. In particular, their radii can be as large 
as the vertical pressure scale height of a PPD. The 3D 
vortices found in our study have a unique property: a 
vortex that grows from a single, local perturbation trig- 
gers a 1^^-generation of vortices that grow at nearby lo- 
cations. The 1^^-generation of vortices grows large and 
triggers a 2^^-generation. The triggering of subsequent 
generations continues ad infinitum. The vortices do not 
advect in the cross-stream direction (radial direction of a 
PPD), but the front dividing the vortex-populated fluid 



FIG. 1: Zombie vortices start near the origin in the x-y plane 
with subsequent generations sweeping outward in x. The vor- 
ticity uoz is red for cyclones; blue for anticyclones; and green 
for Uz = 0. This Couette flow has //iV = 1 and a/iV = -3/4. 
The perturbing vortex at the origin cannot be seen because 
the plane shown here is at z = —0.404. The x-y computa- 
tional domain is \x\ < 4.7124; \y\ < 2.3562, and is larger than 
shown. Movies of Figs. 1 and 2 are online [8]. a) t = 64/A/'. 
b) t = 256/A^. c) t = 576/A^. d) t = 2240/A^. See text for 
details. 



from the unperturbed fluid does. This is shown in Figs. 1 
and 2. Because the vortices grow large and spawn new 
generations of vortices that march across the domain of 
a dead zone, we refer to vortices that self -replicate in this 
manner as zombie vortices. To understand star forma- 
tion, it is important to determine how initial noise cre- 
ates instabilities, but here we focus on the formation and 
self-replication of zombie vortices. 

The simplest flow that is not linearly unstable in which 
zombie vortices occur is a vertically stably-stratified 
Boussinesq fluid in an unbounded plane Couette flow. 
The unperturbed velocity is v = V{x) y with V{x) = ax, 
where a is the uniform shear, x and y are the cross-stream 
and stream-wise coordinates respectively, and quanti- 
ties with "hats" are unit vectors. The span-wise direc- 
tion is which is also the direction of density strati- 
fication and gravity so that the unperturbed density 
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FIG. 2: Zombie vortices sweep outward from the perturbing 
vortex at the origin in the x-z plane shown at 2/ = 0. Anti- 
cyclonic uoz is black and cyclonic is white. This is the same 
flow as in Fig. 1. The computational domain has \z\ < 4.7124 
and is larger than shown. The Rossby number of the initial 
perturbing anticyclone at the origin is Ro = oJz/ f = —0.31. 
a) t = 128/ N. Only the vortex at the origin is present, but 
critical layers with s = and \m\ = 1, 2, and 3 are visible. 
The faint diagonal lines correspond to internal inertia- gravity 
waves with shear, not critical layers, b) t = 480/ N. 1**- 
generation vortices near |x| = 1 and 1/2 have rolled- up from 
critical layers with s = and |m| = 1 and 2, respectively, c) 
t = 1632/Ar. 2'''^-generation \m\ = 1 vortices near |x| = 
and 2 were spawned from the 1^* generation vortices near 
|x| = 1. Another 2^^-generation of \m\ = 1 vortices is near 
|x| 1/2 and 3/2, which were spawned by the 1** generation 
near |x| = 1/2. All 2^^-generation vortices rolled up from 

and 3"^ 



critical layers with \m\ = 1. d) t = 3072/Ar. V\ 2^ 
generation vortices. 



is p{z) = po(l — N'^z/g)i where po is a constant and 
= — g(dp / dz) / p{) is the Brunt- Vaisala frequency, 
which is uniform throughout the domain. The flow is 
on a rotating turntable with angular velocity l^z = //2z, 
where / is the Coriolis parameter. 

Critical layers. In the rotating frame, the dissipation- 
less equations for the velocity v are: 
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where 11 is the pressure head. When Eqs. ([I]) - (|3]) are 
linearized about v = V{x) y and the eigenmodes are 
proportional to ^"^{.Kv+^zz-st) ^ where we assume periodic 



boundary conditions in z and y. When the density p is 
stably-stratified (i.e., with real 7^ 0) or constant (A^ = 
0) this plane Couette flow is neutrally stable for all Vt (i.e., 
the real part of s is zero and the eigenmodes neither grow 
nor decay). The eigen-equation for the stream functions 
of the eigenmodes of Eqs. ([1]) - ([3]) is a generalization 
of Rayleigh's equation [9] and is a 2^^-order o.d.e. The 
coefficient of the highest-derivative term is 



[V{x)-S/ky]{[V{x)-S/ky] 



[N/kyf]. (4) 



It is well known that the eigenmodes of an o.d.e. are sin- 
gular at locations where the coefficient of the highest- 
derivative term of the eigen-equation is equal to zero. 
At X*, the eigenmode has a critical layer [9[. For flu- 
ids with constant density (A' = 0), the critical layers of 
uni-directional shearing flows have been well-studied and 
have y(x*) equal to the phase speed s/ky of the eigen- 
mode. We refer to these as harotropic critical layers; they 
are not of interest to us because laboratory experiments 
and numerical calculations show that these critical lay- 
ers are not easily excited and do not produce vortices. 
Eigenmodes with barotropic critical layers have singular- 
ities in the stream-wise components of their velocities, 
but not other components. 

For A^ 7^ 0, expression (j4]) shows that there are eigen- 
modes with barotropic critical layers at V{x*) — s/ky, but 
we found that they remain difficult to excite and never 
form vortices. However, there is a new class of eigen- 
modes with critical layers with V{x^) — s/ky^N/ky = 0. 
We refer to the latter as baroclinic critical layers. These 
have not been explored before. Eigenmodes with baro- 
clinic critical layers are neutrally stable (i.e., s is real). 
These eigenmodes are easily excited by small distur- 
bances, always produce large-amplitude vortex layers, 
and often produce zombie vortices because their singular- 
ities are in the vertical z components of their velocities. 
To see why, note that the z-component of the curl of 
Eq. ([I]) gives 

du./dt = -(v • + (u; • V)^, + (/ + a){dv,/dz), (5) 

where u; is the relative vorticity with respect to the equi- 
librium flow, i.e., = V X (v — V{x)y). Equation (|5]) 
shows that the generalized Coriolis term {f -\-a){dvz/dz) 
is a source for cj^. Within the baroclinic critical layer, 
the z-component of the velocity is nearly anti-symmetric 
about the plane with x = x*; on one side of the layer 
Vz ^ 00^ and on the other side Vz —00; thus, the 
CorioHs term in Eq. ^ creates a large-magnitude vortex 
layer centered at x* made of dipolar segments with one 
side having cyclonic vorticity {oJzf > 0) and the other 
with anticyclonic vorticity {oJzf < 0) (cf., Fig. 1(a)). 

We consider only the case of anticyclonic shear with 
a < and / > 0. Note that PPDs have anticyclonic shear 
and are periodic in their stream-wise, or azimuthal, di- 
rection like the Couette flow studied here. The equations 
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of motion (1) - (3) and boundary conditions are invariant 
under translations in y and and also under translation 
in X by (5 when accompanied by a stream-wise boost in 
velocity of (j5. The latter symmetry is known as the shift 
and boost symmetry^ c.f., [iqI^IiiI, and is exploited when 
using shearing sheet boundary conditions of PPDs [HI . 
From this point on, we use non-dimensional units with 
time in units of 1 /N and length in units of | {LN) / (27rcr) | , 
where L is the periodicity length in the stream-wise di- 
rection, so ky in expression (|4]) is 27rm/I/, where m is 
an integer. Baroclinic critical layers have ky ^ 0, and 
expression (|4]) shows that in dimensionless units are at: 

x* = -{s±l)/m. (6) 

Due to the shift-and-boost symmetry, the origin of the x- 
axis is not unique, so Eq. ([6j) has the following meaning: 
X* is the cross-stream distance between a finite- amplitude 
perturbation and the location of the baroclinic critical 
layer that it excites. 

Single vortex perturbation. Our numerical solutions 
of Eqs. (1) - (3) show that many types of perturbations 
create zombie vortices. For verification, our solutions 
were computed with two independent codes using the 
methodology in [l2|. One code treated the cross-stream 
boundaries as periodic with energy damping, and the 
other used the shearing sheet approximation. In the ex- 
ample in Figs. 1 and 2, the initial perturbation is a small- 
volume anticyclone with a Rossby number Ro = ujz/ f oi 
—0.31 at the origin superposed on the unperturbed fiow 
V{x) and p{z). Figure 1 shows uOz in an x-y plane. The 
perturbing vortex is steady, so it excites critical layers 
with temporal frequencies 5 = 0. Thus, Eq. (|6]) shows 
that the critical layers are at \x*\ = l/\m\ and that the 
vortex at the origin has no critical layers with x > 1. At 
early times. Fig. 1(a) shows that vortex layers form at 
the critical layers: at each critical layer with wave number 
|m|, uoz appears at x = l/\m\ as |m| segments of dipolar 
stripes aligned in the stream-wise y direction. A Fourier 
analysis shows that the stripes have s = 0. We have 
shown previously ll3|,ll^ that in shear fiows with /a < 0, 
cyclonic vortex layers aligned in the stream- wise direction 
are stable, whereas anticyclonic layers are unstable, roll- 
up into discrete anticyclones, and merge to form one large 
anticyclone. This behavior is seen in Fig. 1(b). The an- 
ticyclonic vorticity at x = 1/3 has already rolled up and 
merged into a single large anticyclone (near y = 1.5). 
The anticyclonic vorticity at x = 1 /2 has rolled up into a 
large anticyclone near y = —0.5. In contrast, the cyclonic 
(Jz near x = 1/2 has formed a continuous, but meander- 
ing, filament. At later times (Fig. 1(c)) the anticyclones 
near x = 1/3 (and near y = 2) and near x = 1/2 (and 
near y = —1) have become larger. Moreover, Figs. 1(c) 
and 1(d) show critical layers and vortices at x > 1, which 
cannot be created by the perturbation at the origin. The 
layers at x > 1 are due to the self-replication of the 1^^- 
generation vortices at \x\ < 1. A vortex at any location 



will excite critical layers in a manner exactly like the 
original perturbing vortex at the origin due to the shift- 
and-boost symmetry (and will have s = when viewed 
in the frame moving with the perturbing vortex). Fig- 
ure 1(c), shows 2^'^-generation critical layers at x = 4/3, 
3/2, 2, and 2/3 ah with |m| = 1 and excited by the 1^^- 
generation vortices at x = 1/3, 1/2, 1, and —1/3, respec- 
tively. Figure 1(d), shows 3^^-generation critical layers 
at 2 < X < 3, and 4^^-generation critical layers becoming 
excited at x > 3. At late times the vortices from m = 1 
critical layers dominate the fiow. (See Fig. 2(d).) At very 
late times, the vortices have widths in their cross-stream 
directions of order unity. (See below.) 

Figure 2 shows the fiow in Fig. 1 viewed in the x-z 
plane and illustrates our main result: at late times the do- 
main in the x-z plane fills with anticyclones. Because the 
fiow is homogeneous with uniform a and A/", the vortices 
form a regular lattice despite the fiow's turbulence. As 
time progresses in Fig. 2, the vortex population spreads 
out from the perturbing vortex at the origin. At early 
times (Fig. 2(a)) the fiow has 1^^-generation critical lay- 
ers, with |m| = 1, 2, and 3 being the most visible. In this 
first generation, and all subsequent generations, a per- 
turbing vortex creates four vortices at its |m| = 1 critical 
layers. The new vortices form at locations in x that are 
:^lx distant from the perturbing vortex and at locations 
in z that are i/^ distant from the perturbing vortex. The 
distance in x between a perturbing vortex and the center- 
line of the |m| = 1 critical layer it excites (and the center- 
line of the dipolar vortex segment it creates) is unity, but 
the distance to the anticyclonic piece of the dipolar vor- 
tex is Ix , which is slightly greater than unity. At the time 
of Fig. 2(b) the 2"^^-generation m = 1 critical layers cre- 
ated by the 1^^-generation vortices with |m| = 1, 2, and 
3 are all visible near \x\ = 2lx, (3/2) Z^^, and (4/3) re- 
spectively. At later times (Fig. 2(d)), the \m\ = 1 vortices 
that were descended from the l^*-generation \m\ = 1 vor- 
tices dominate the fiow. These dominant vortices form 
a lattice and are located at [x = 2n Ix^ z = 2j Iz] and at 
[x = {2n + l)lx-,z = (2j + 1)1 z\i for all integers n and j. 

The characteristic \Ro\ of late-time zombie vortices is 
of order that of the initial perturbation. After a vortex 
forms, its \Ro\ intensifies to its approximate peak value 
within a few of its turn- around times, and it remains 
near that value indefinitely. Based on several numerical 
experiments, it appears that the late-time values of \Ro\ 
depend on the parameters, f and a rather than on 
the properties of the initial perturbation. To examine the 
energy of the vortices and discover its source, we decom- 
posed the fiow's energy into two orthogonal parts: (1) 
the zonal component consisting of the kinetic energy of 
the stream-wise velocity component with Fourier modes 
kz = ky = (i.e., the background shearing fiow); and (2) 
the non-zonal component consisting of everything else, 
including the potential energy g J z{p — p) (d volume). If 
the initial fiow were the unperturbed fiow V{x)y and p. 
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then the initial energy would be all zonal. In the flow in 
Figs. 1 and 2, there is a small initial non-zonal compo- 
nent of the energy due to the initial vortex at the origin. 
At late times, the non-zonal energy component repre- 
sents the energy of the zombie vortices (and their tur- 
bulence and waves). After initial transients with faster 
than exponential growth, the non-zonal energy increases 
exponentially from just after the time of Fig. 1(a) to the 
time in Fig. 2(d) with an e- folding time of ^ 1000. The 
non-zonal energy component in Fig. 2(d) is more than 
400 times larger than its initial value, and that energy 
is supplied by the zonal component. The total energy 
is not conserved either with shearing sheet boundaries 
or with energy-damping periodic boundaries. In the lat- 
ter, the vortices primarily lose energy by radiating in- 
ternal, shearing inertia-gravity waves to the boundaries. 
If the self-replication is self-similar, we would expect the 
perimeter of the front between the vortex-dominated flow 
and unperturbed flow in each x-z plane to grow as t and 
the number of vortices to increase as t^, which is consis- 
tent with our numerical experiments. After transients, 
the energy of the vortices, i.e., the non-zonal energy com- 
ponent, grows exponentially in time because the volume 
of each vortex grows exponentially in our calculations. 
That exponential growth must eventually stop because 
the cross-stream widths of the vortices cannot become 
larger than unity without them overlapping. 

Discussion. Neutrally-stable, shearing, stably- 
stratified Couette flow can be finite- amplitude unstable. 
In the example here, baroclinic critical layers are excited 
by a small vortex, but our calculations show that a 
variety of small- volume, small-energy perturbations 
cause critical layers to grow and roll-up into large- 
volume, large-energy vortices. In general, this instability 
self-replicates with each new vortex exciting new layers 
that roll-up until the domain fills with compact 3D 
(i.e., not Taylor columns) vortices. The robustness of 
zombie vortices is self-evident from the fact that they 
are embedded in turbulent flows and survive indefinitely. 
They survive by drawing energy from the background 
shear flow at a rate that is faster than their dissipation 
rate. For constant N and a, the unperturbed flow is 
homogeneous, and vortex self-replication is self-similar 
with zombie vortices forming a regular lattice. The 
regularity of the lattice allows for reinforcement: each 
vortex re-excites four other vortices in the lattice, and 
each vortex in the lattice is continually re-excited by 
four other vortices. Zombie vortices occur frequently 
in our numerical simulations of Boussinesq and com- 



pressible fluids, so they pose a paradox: if they are so 
common, why have they not been reported earlier? We 
believe there are three reasons: (1) instabilities have 
not been systematically sought in stratifled Couette 
flows (2) with few exceptions [15j, stability studies 
of ideal gases in PPDs were carried out with no initial 
vertical stratiflcation jl, and (3) the necessary spatial 
resolution to compute critical layers is lacking in many 
calculations. Zombie vortices occur in our calculations 
of the dead zones of protoplanetary disks [16], which 
suggest that they may have an important role in star and 
planet formation. In addition, zombie vortices should 
be observable in laboratory circular Couette flows with 
stratifled salt water for parameter values where the flow 
is linearly stable and where the angular momentum of 
the flow increases with increasing radius. 
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